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Abstract 


In this paper, we focus on investigating a post-processing technique de- 
signed for one-dimensional singularly perturbed parabolic convection-diffusion 
problems that demonstrate a regular boundary layer. We use a back- 
ward Euler numerical approach for time derivatives with uniform mesh 
in the temporal direction, and a simple upwind scheme is used for spa- 
tial derivatives with modified graded mesh in the spatial direction. In 
this study, we demonstrate the effectiveness of the Richardson extrapola- 


tion technique in enhancing the c-uniform accuracy of simple upwinding 
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within the discrete supremum norm, as evidenced by an improvement from 
O(N~! In(1/e) + A@) to O(N~? In?(1/e) + A@?). Furthermore, to validate 
the theoretical findings, computational experiments are conducted for two 


test examples by applying the proposed technique. 
AMS subject classifications (2020): 65M06, 65M12, 65M15, 65M22. 


Keywords: Singularly perturbed parabolic problem; Regular boundary layer; 
Upwind scheme; Richardson extrapolation; Modified graded mesh; Uniform 


convergence. 


1 Introduction 


We consider the following one-dimensional singularly perturbed parabolic 
convection-diffusion initial-boundary value problem (IBVP) posed on the do- 
main N = A, x Ag, where A, = (0,1), Ag = (0,7): 


Ley(r8) = (AHEM PME) 5g (ry MEO + walryuir6)) 


= f(r,9), (r,0) ER, 
y(r,0)=yo(r), rey, 
y(0,0) = 0, 0 € No, 
y(1, 0) = 0, 6 € No, 


(1) 
where 0 < € < 1 is a small parameter and the coefficients V; and Wz are 


sufficiently smooth function such that 
Wi(r)>2a>0, Wa(r)>B>0 on Ay. (2) 


Under sufficient smoothness and compatibility conditions imposed on the 
functions yo and f, the parabolic IBVP (1)—(2), in general, admits a unique 
solution y(r,@), which exhibits a regular boundary layer at r = 1. This 
type of problem encompasses the linearized Navier-Stokes equation, which 
emerges in the modeling of convection-dominated flow issues in fluid dynam- 


ics, particularly when dealing with large Reynolds numbers. In the presence 
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of a boundary layer, conventional numerical methods, such as standard fi- 
nite difference or finite element schemes, when applied on uniform meshes, 
fail to produce accurate numerical solutions as the singular perturbation pa- 
rameter € approaches zero. This limitation serves as a driving force to devise 
é-uniformly convergent numerical methods. Among these approaches, the fit- 
ted mesh method stands out as a satisfying and widely embraced technique 
for surmounting numerical challenges. By employing an especially designed 
layer-adapted mesh, this method successfully overcomes the shortcomings 


encountered with traditional approaches. 


Over the past few years, numerous authors have made significant strides 
in the development of uniformly convergent numerical methods on Shishkin 
meshes, specifically tailored for singularly perturbed parabolic convection- 
diffusion problems. These methods address the challenges posed by such 
complex scenarios and have shown great promise in achieving convergence, 
mainly Cai and Liu [1], Clavero, Jorge, and Lisbona [3], O’Riordan, Pick- 
ett, and Shishkin [28] in the presence of regular boundary layers. However, 
despite these advancements, it is important to note that all of these meth- 
ods exhibit first-order accuracy in both the spatial and temporal variables. 
Hemker [9] employed a defect correction technique to enhance the accuracy 
of temporal variable computations for parabolic singularly perturbed prob- 
lems, specifically those lacking the convective term. Additionally, in [10], the 
same approach was applied to parabolic convection-diffusion problems. This 
technique has proved to be effective in refining the precision of the results in 


both cases. 


Woldaregay and Duressa [35] considered the study of singularly perturbed 
differential-difference equations having delay and advance in the reaction 
term. Debela, Woldaregay, and Duressa [6] studied the singularly perturbed 
differential equations of convection diffusion type with nonlocal boundary 
conditions. Mekonnen and Duressa [15] considered the study of numeri- 
cal treatment of two parametric singularly perturbed parabolic convection- 
diffusion problems. The scheme is developed through the Crank—Nicholson 
discretization method in the temporal dimension followed by fitting the B- 


spline collocation method in the spatial direction. Also, Hailu and Duressa [8] 
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focused on the study of singularly perturbed parabolic convection-diffusion 


equations with integral boundary conditions and a large negative shift. 


Furthermore, Natesan and Deb [5] have recently introduced two novel 
numerical schemes that exhibit uniform convergence and higher-order accu- 
racy in time for singularly perturbed parabolic reaction-diffusion problems. 
These schemes are specifically designed to handle scenarios involving only 
parabolic layers. In their recent work, Mukherjee and Natesan [17] proposed a 
novel hybrid numerical scheme for singularly perturbed parabolic convection- 
diffusion problems. This innovative approach demonstrates uniform conver- 
gence, achieving an impressive order of convergence one in time and nearly 
two in space. The results of their study pave the way for more efficient and 
accurate solutions to such challenging problems. Several researchers, includ- 
ing Clavero, Gracia, and Jorge [2], Kopteva [14], and Shishkin [30], have made 
significant contributions to the development of uniformly convergent methods 
of order two in both variables for singularly perturbed parabolic convection- 
diffusion problems. In their study, Mukherjee and Natesan [18] employed 
the Richardson extrapolation technique to enhance the ¢-uniform accuracy 
of the simple upwinding method. This improvement was observed in the 
discrete supremum norm and extended the accuracy from O(N~'In N 4+ At) 
to O(N~* In? N + At?) on the Shishkin mesh. The focus of their investiga- 
tion was one-dimensional singularly perturbed parabolic convection-diffusion 


problems that exhibit a regular boundary layer. 


Maneesh and Natesan [33] presented an advanced numerical approach 
designed to handle singularly perturbed systems of parabolic convection- 
diffusion problems, which feature overlapping exponential boundary layers. 
The proposed method achieves uniform convergence of higher-order accu- 
racy. The numerical scheme integrates the backward-Euler method for the 
time derivative on a uniform mesh, along with the classical upwind scheme 
for spatial derivatives on a piecewise-uniform Shishkin mesh. This combina- 
tion yields near first-order convergence in both space and time. To further 
enhance precision, the authors employed the Richardson extrapolation tech- 
nique, elevating the accuracy of the scheme to second-order, maintaining 


uniform convergence in both time and space. 
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Negero [20, 25] considered the study of an exponentially fitted numeri- 
cal scheme and analyzed it for solving singularly perturbed two-parameter 
parabolic problems with large temporal lag. Negero [22] studied the numer- 
ical approximation for two-parameter singularly perturbed parabolic partial 
differential equations with time delay. In [23], the author considered the 
class of time-delayed, singularly perturbed parabolic reaction-diffusion prob- 
lems. Negero [21, 24, 26] studied time delay problems and parameter-uniform 
robust schemes for singularly perturbed parabolic convection-diffusion prob- 
lems with large time-lag. 


In this current study, our objective is to employ and examine a straight- 


forward post-processing technique on the basic upwinding solution. The 
primary goal is to attain a convergence order greater than one concern- 
ing both the spatial and temporal variables. Richardson extrapolation, a 
renowned post-processing technique, offers a superior approximation to the 
exact solution achieved by averaging the numerical solutions computed on 
two embedded meshes. This method has been extensively investigated in 
the literature [11, 31] to enhance the accuracy of numerical solutions for sin- 
gularly perturbed elliptic reaction-diffusion equations. However, the crux of 
their analysis relies on the direct expansion of the upwinding solution, which 
proves to be intricate. Shishkin and Shishkina [32] have consistently adopted 
this approach for quasilinear parabolic convection-diffusion problems. How- 
ever, Natividad and Stynes [19] presented a comparatively simpler and dis- 
tinctive analysis for applying Richardson extrapolation to one-dimensional 
singularly perturbed convection-diffusion boundary-value problems (BVPs). 
More recently, Deb and Natesan [4] conducted an analysis of the Richardson 
extrapolation technique for solving singularly perturbed coupled systems of 


convection-diffusion BVPs. 


This paper represents the inaugural analysis of Richardson extrapolation 
on the nonuniform graded mesh applied to singularly perturbed parabolic 
convection-diffusion IBVPs through error decomposition after extrapolation. 
Initially, we address the IBVP (1)—(2) on a nonuniform graded mesh using the 
classical implicit upwind finite difference scheme. Subsequently, we demon- 
strate that the implicit upwind scheme achieves e-uniform convergence with 


an almost first-order accuracy in the discrete supremum norm. Next, we 
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incorporate the Richardson extrapolation technique to enhance the nearly 
first-order convergence achieved by the simple upwinding method, leading 
to an almost second-order convergence. During the error analysis, we derive 
separate estimates for the smooth component and the layer component of the 


error. Consequently, we attain the necessary ¢-uniform convergence result. 


The subsequent sections of this paper are organized as follows: In sec- 
tion 2, an apriori bound on the analytical solution is presented, and stronger 
bounds on the derivatives of the solution are derived through decomposition. 
Section 3 introduces the modified graded mesh and the classical implicit up- 
wind scheme used for discretizing the continuous problem. Here, we establish 
the e-uniform convergence result for the implicit upwind scheme and provide 
some technical results to be utilized later in this paper. Moving on, section 4 
introduces the Richardson extrapolation technique, and the main theoretical 
result, demonstrating that the extrapolated solution of the upwind scheme 
achieves ¢-uniform convergence with almost second-order accuracy to the ex- 
act solution of the continuous problem, is proven. In section 5, numerical 
experiments are conducted on two test examples to validate the theoretical 


results. Finally, section 6 summarizes the main conclusions of this study. 


In this paper, we use the notation C to represent a universal positive 
constant that remains unaffected by changes in the perturbation parameter 
€, the values of N and M (representing the number of mesh intervals in the 
spatial and temporal directions, respectively), as well as the mesh points. In 
the analysis, we use the standard supremum norm || - ||..,p, which is defined 
by 


I|Zll0,p = sup |2(x)], 
KED 


for a function z defined on some domain D. 
Before we analyze the problem, some of the compatibility conditions are 
necessary. Therefore, the following compatibility conditions at the corners 


for functions and its first-order derivatives are assumed to satisfy, 
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Then (1) has a unique solution in the parabolic Holder space O?+%1+¢/2(X); 
see [29, 7]. Moreover, the fulfillment of second-order corner compatibility 
conditions requires the Holder space to be C4+%?+¢/?(&). These conditions 
can be written down explicitly in terms of the data of the problem in the 


following way: differentiating (1) with respect to 0, we get 


fo = yoo + Leyo + Vooy = yoo + Le(f — Ley) + Vaoy. 


Hence, recalling (1) and (3), the second-order corner compatibility condition 
is 

L-(Leyo) = Lef — fo (4) 
at the corners (0,0) and (1,0). Under these hypotheses, the solution y of (1) 


has an exponential layer along the boundary r = 1 of X. 


2 Bounds on the solution and its derivatives 


In this section, we introduce the conventional a-priori bound for the analytical 
solution and further establish enhanced bounds on the derivatives of the 
solution for problems (1)—(2). This is achieved through a decomposition of 
the solution into smooth and layer components. The utilization of these 
bounds becomes essential in the subsequent section as we aim to prove the 
e-uniform error estimate. We first present some initial reports. Let G = 8/X. 
The differential operator D- then fulfills the following minimum principle on 
XN, and [29, Theorem 2.2] provides the evidence for this. 


Lemma 1. [Minimum principle] Assume that a function z € C°(&) N C?(X) 
satisfies z(r,0) > 0,(r,0) € G and Lez(r,@) > 0,(r,0) € &. Then we have 
2(r,0) > 0, for all (r,0) ER. 


The following is a direct result of the aforementioned minimum principle: 


é-uniform bound of the solution of the problem (1)—(2). 


Lemma 2. The solution y of the IBVP (1)—(2) satisfies 


ly(r,A)) SC, (7,0) ER. 
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Proof. The proof follows from Lemma 2.3 of Roos, Stynes, and Tobiska [29]. 


Subsequently, we present a priori bounds concerning the derivatives of the 
analytical solution y to problem (1)—(2) in the forthcoming theorem, which 


will be utilized in Theorem 2. 


Theorem 1. For all nonnegative integers p,q satisfying 0 < p+q < 5, the 
exact solution y of the IBVP (1)-(2) satisfies the estimate 
| oPptdy 


arPaed (r, 0) < o(1 +e Pexp(—a(1— nie), (r,0) ER. (5) 


Proof. In [27], this conclusion was verified for 0 < p+q < 2,. Under necessary 


compatibility conditions and sufficient smoothness on the data, the proof of 


the estimate (5) for higher values of p,q follows similarly from [2]. 


Let us now decompose the solution y of the IBVP (1)-(2) into the form 
of y=1+™m, where | and m represent the smooth component and the layer 
component, respectively. We further break down the smooth component into 


the sum 
4 


i=) le 


j=0 


where the functions /;, 7 = 0,1,2,3, are solutions of the following first-order 


problems: 
dlo , Slo 7, 2 
a6 t Wy Ap T Wolo f in X, 
1o(0,8) = y(0,8), 8€(0,T], lo(r,0) = y(r,0), re Ar, 6) 
Al; al; — Ohya , 
ya + barr + Wal; —_ Or2 m X, 


[,(0,0)=0, @€(0,T7], J;(r,0)=0, reA,, j=1,2,3, 


and the lastly, the function l, satisfies 


ls 
Lel4 = Dre in X, (7) 


14(0,0) =14(1,0) =0, 0€(0,T7], lu(r,0)=0, re A,. 


Hence, the smooth component / satisfies the following IBVP: 
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Ll=f in X, 
(0,6) =0, 1(1,8) = Y5_9€7j(1,0), 8 € (0,71, (8) 
i(r,0) =y(r,0), re Ay, 

and therefore, the layer component m must satisfy 
L-m=0 in X, 
m(0,@) =0, m(1,@) = y(1,@) — 11,6), @ € (0,7), (9) 
m(r,0) =0, reéA,. 

Theorem 2. For any nonnegative integers p and q that fulfill the condition 


0 < p+q < 5, the smooth component / and the layer component m, as 


defined in (8) and (9), respectively, adhere to the subsequent bounds, 


SOs ey, 


00,8 


OP+4] 
| OrP 004 


OPtIm 

Serooe 8) < Ce Pexp(—a(l—r)/e), (r, A) ER. 

Proof. In the initial step, we will establish more stringent bounds for the 
smooth component 1, as defined in (8), and its derivatives. The functions 
lj,j = 0,1,2,3, serve as solutions to the problems outlined in (6). Impor- 
tantly, these functions remain unaffected by the parameter ¢. As a conse- 
quence, their derivatives exhibit the property of ¢-uniformly bounded behav- 


ior. Hence 


<C, j=0,1,2,3, for O0<pt+q<5. (10) 


arta, 
| OrP 064 


o0,N 
Similarly, the function l, corresponds to the solution of a problem akin to (1). 
By applying estimate(5) in a manner analogous to l4, we obtain the bounds 


<C(l+e ?exp(—a(1—r)/e)), for O<p+q<5. (11) 


OP+4] 4 
| OrP 004 


o00,N 


Therefore, by combining the estimates (10) and (11), for 0 < p+q< 5, we 


can establish the necessary bounds on the smooth component / as follows: 
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Orta, 
OrP 064 


< J 


& 
00,8 j=0 


| Orta] 


Drage <C(l +e"). 


o0o0,N 


Alternatively, we can utilize the minimum principle (Lemma 1) on X with 


the barrier function 
W(r, 0) = Cexp(—a(l—r)/e), (7,0) EX. 


By selecting a sufficiently large C, we achieve the necessary bound on m. 


Moreover, the bounds on the derivatives of m can be deduced, following the 


arguments presented in [16, 28], thus concluding the proof. 


3 Numerical discretization 


Within this section, a well-suited mesh is presented for discretizing the 
domain, ensuring the attainment of an e-uniformly convergent difference 
scheme. Additionally, a detailed exposition of the employed difference scheme 


for discretizing the problem (1)—(2) is provided. 


3.1 Temporal Discretization 


To establish the convergence of (1)—(2) at each instance, we use the uniform 


time grid 
wh = {0, =nA0, n=0, 1, ..., M, AO=T/M}, 


Here M represents the grid points. 


3.2 Spatial discretization 


We generate a modified graded mesh, AN in the interval [0,1] and order to 


resolve the boundary layer at r = 1 as follows: 


pg =1=8n-1 for J, (12) 
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where ju is defined as the following form: 


Vo = 0, 
0; =2e4, 1<j <4, 
J N = 2 (13) 
Oj =0;(1+ ph), F<j<N-2, 
Un =1, 
where the parameter h satisfies the following nonlinear equation: 
In(1/e) = (N/2)In(1 + ph). (14) 


The above section of the parameter h ensures that there are N/2 grid points 
in the interval [0,1 — ¢] which are distributed gradely in the interval [0,1 — 
e]. Numerical verification stimulate that the interval (Jy_1,1) is not too 
small in comparison with the previous one (Wy_—2,Vn-1). In the subinterval 
[1 — ¢,1], we distribute N/2 points with uniform step length 2¢/N, while 
in the subinterval [0,1 — ¢], we first find h for some fix N, by means of 
the nonlinear equation (14), and corresponding to that h we distribute N/2 
points in the interval [0, 1—¢]. The mesh length is denoted by hj = 0; —V;-1, 
for 7 =1, 2, ..., N. 


Remark 1. The mesh size in piecewise uniform and the modified graded 


region is given by 


2e/N for j=1, 2, ..., N/2, 
phd;—-1 for j = N/2+1, N/2 +2, aay dN 


Lemma 3. The mesh defined in (13) satisfies the following estimates: 


for j=1, 2, ..., N/2, 
|rj+1 — hyl < 
Ch for j=N/24+1, N/24+2, ..., N. 
Proof. Initially, we consider 7 = 1, 2, ..., N/2. As the mesh is uniform in 


this portion, so there is nothing to prove. 
For j7 = N/2+1, N/2+2, ..., N, we have 


|njz1 — hy] = |phd; — phd;-1| 
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Figure 1: Modified graded mesh layer at r = 1 for values of N and ¢ (a) N = 32,e= 1071 
(b) N = 128,e= 1072. 


= ph|d; — 0;-1| 
= prhojz1 
<Ch. 


Here, we have taken 0 < p,h <1. 


Lemma 4. For the modified graded mesh defined in (13), the parameter h 


satisfies the following bound: 
h< CN In(1/e). 


Proof. Let Ky, be the number of points J; in the partition (13) such that 
0; <e, forj =1, 2, ..., N/2. Clearly kK, < C/h and Ko are the number 
of points in the partition (13) such that J; > e. Let Uy/241 be the smallest 
point such that ¥; > ¢. We have to estimate the bound for Ky. Assuming 
ph <1, we have 


a os Oj41 
iS So i= Oi —9,)7 f dd 
N/241 N/241 05 
ad Oj+1 
-> (i f dd 
N/241 95 
iN. Oj+1 
= > (ona fa 
N/241 9; 
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< Y @/phdjay? fav, 
N/241 5 


because Jj+1 < 20,;. For any J € [J;,V41], we have 


N itl | 
Ko< > 2(ph)-* f = 
N/2+1 95 


Recalling N = Kk, + Ko, we have 
N < C/ph+2(ph)~* In(1/e) 
N <1/h(Cp + 2(p)~* In(1/e)) 
N <1/h(C In(1/e)). 


Finally, we get 
h<CN'In(1/e), 


where N is the number of grid points in the r-direction. 


3.3 The classical implicit upwind scheme 


Prior to explaining the scheme, we establish necessary operators for a given 
mesh function z(r;,9,), in space and time. First, we define the forward dif- 
ference operator D;*, the backward difference operator D7, and the central 
difference operator D® in the spatial domain. These operators allow us to 
approximate derivatives with respect to the spatial variable, enabling us to 
capture the behavior of the function at different points. Additionally, we 
introduce the backward difference operator Dy, in the temporal domain rep- 


resented by 
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Dt 2lr;, On) = 2(75341;9n) _ 2(rj, On) 


j Ajsa ; 
De z(1j,) = “iz bn) os nm) 
z 15 
D8 z(r3, On) = ae ar 5-1 On) ( ) 
J 
i 2z(rj,On — 2(75, On— 
Dg 2(7j,9n) = (73 On) = j 1) 


respectively, and define the second-order finite difference operator 6? in space 


by 


She Pipe ee ee Pe) 
r 239° Rp x ; 


hj 
where h; =hy+hj+1 j=1,...,N—-—land hj =; —T5j-1; FS Tyg lV 
Let xVM = nM C8 and let TY! = ae \ RY. For the discretiza- 


tion of the continuous problem (1), the following classical implicit upwind 


; : : <N,M 
finite difference scheme is used on the mesh XN, : 


LNMYNAO = (Dg — 662 + UW D> + W2)Y"? = f in NVM 


(16) 
YA =y, on GYM 
’ eS - 
Demonstrating the discrete minimum principle of the finite difference op- 
erator fia *M establishes its well-known property, ultimately leading to the 


e-uniform stability of the difference operator LN. 


Lemma 5. [Discrete Minimum Principle] Assume that a mesh function 

Usatisfies U >00n GY™. Then LN™“U > 0 in XY implies that U > 0 at 
<N,M 

each point of RN.” . 


Proof. Let s and | be indices such that U! = min; ,for U? € RY Assume 
J, 


that U! < 0. It is easy to see that (s,l) € {1,2,...,N} x {1,2,...,M}, 
because otherwise U! > 0. It follows that Ul,, —U! > 0 and U!_, —U! >0. 
Thus LN“U! < 0, which is a contradiction. Therefore U! > 0. The indices 


=) 
s and | being arbitrary, we obtain Ue 20 ink;. 


The subsequent theorem asserts the ¢-uniform convergence of the numer- 
<N,M 
ical scheme (16) on the modified graded mesh XN.” , demonstrating nearly 


first-order accuracy. 
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Theorem 3. [Error due to up-winding] Consider the solution y to the con- 
tinuous problem defined by (1)-(2), and let Y“:4® represent the solution 
corresponding to the discrete problem outlined in (16). We can then exam- 
ine the error linked to the discrete solution Y’4? at the time instance 6,, 
that satisfies 


y(75,9n) — YN? (r;,On)| < C(N7*In(1/e)+ AO), forl<j<N-1. 


Proof. We obtained this result with the help of article Mukherjee and Nate- 


san; see [18, Appendix A]. 


The primary objective of this article is to enhance the accuracy of the 
discrete solution YN“? for the problem (16) through a post-processing tech- 
nique. The aim is to achieve an ¢-uniform order of accuracy greater than one, 
considering both the spatial and temporal variables for the IBVP (1)-(2). 
To accomplish this, we will employ the Richardson extrapolation technique. 
Prior to introducing this technique, we present some essential technical lem- 


mas that will be utilized in section 4. 


Lemma 6. On AY = {r;}’, define the mesh function 


’ C1 
If Re ( -) (17) 
for some constant C1. Moreover, for N/2+1< j < N —1, we have 
LE Re > Coe" RB, (18) 


for some constant C. 


ie 
Proof. Firstly, Rj — Rj-1 = ack ae we have 
€ 


2a a 
LNMR. — R;—R,_ W;—R;_1+U;R; 
E J aa ju) + Te ae 


wes oe W; — 2ahj | 
€ h; + hg4i 
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Ca 


ey 1<j<N-1. 
~ et+ah; for j 


Therefore, (17) is proved and (18) is an easy consequence of it, since h/e < 
A/a. 


Lemma 7. For the modified graded mesh rd = {r;}q’, the following in- 


equalities hold true: 
(i) 
4, hte =a 
exp(—a(1—1;)/e) < IL (+ am) , forallj. (19) 


(ii) There exists a constant C' such that 


ee ’ 
II (1 4 =~) <CN-*0-3/§) for N/2<j7<N-1. (20) 
€ 


Proof. The proof of (2) and (#) follows from [834]. 


4 Extrapolation of Y"’*? 


Within this section, we introduce the extrapolation technique and concur- 
rently present a thorough error analysis. This analysis involves the decom- 
position of the solution pertaining to the discrete problem outlined in (16). 
Ultimately, we establish an ¢-uniform error estimate that is closely linked 


with the extrapolated solution. 


4.1 Extrapolation technique 


To improve the accuracy of the numerical solution Y%’4? by extrapolation, 


2N,2M Z2N,e 


we shall solve the discrete problem (16) on the fine mesh Ne =A xX 


W2™ , where A, = (0,1) with 2N mesh interval in the spatial direction and 
2M mesh interval in the 6-direction, where aaa is a Modified graded mesh 
having the same transition point (1—¢) as 7a * and is obtained by bisecting 


each mesh interval of ree Clearly, RM ={(r7, Ganj be Vea ae = {(¥;,4n)} 
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and the following are the corresponding mesh widths in eNeM . 


H/2, for #; en (0,1 - 4], where H = 2¢/N, 


73—Tj-1=% h/2, for *,€ ia M{l-e,1], where h = phdj-1, 
6,, —On_1 = AO /2 for On € we". 
Let Y?"-49/2 be the solution of the discrete problem (16) on the mesh 
oo , Now, it follows from Theorem 3 that 


YN 49 (rs On) — y(ts,On) = C3 (w- In(1/e) + 49) + Ry,ao(%j, On) 


= C3 (1v-*(ae/22) + 40) + Ry,ao(r;,9n), 


(73,n) ER, (21) 


where C’3 is some fixed constant and the remainder term Ry, Ag is 
O (Ww In(1/e)+ 40) . keeping in the mind that Y2"49/? is obtaining using 


the same transition point 1 — ¢, we have 


Y2NAF/2 (5, On) — y(Fj, On) = C3 (ce) (ae/22) 4 (A6/2)) + Row,ao/2(7j,n), 


ao % S2N,2M 
(Tj, On) € a. ; 


(22) 


where the remainder term R2y,Ag/2 is o( x-*n(a/e) + A6}. Now, the 
elimination of the O(N~') and O(A@) terms from (21) and (22) lead to the 


following approximation: 


y(r;; ce) = (ann, On) _ beac ee (ate on)) 


= —2Row,a/2(rj,On) + Rv,ao(rj, On) 
= o(w In(1/e) + 48), (r;,0n,) ERO 


Hence, we will employ the subsequent straightforward extrapolation formula, 


YN? (5, On) = 2V2NA9/2 (5, n) —YNA%(05, On), (77,9n) € Re, (23) 


e 
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which is expected to yield an approximation of y(r;,0,) is more accurate 


than either YN-49(r;,0,) or Y2N-49/2 (rj, On). 


4.2 Solution decomposition 


To obtain the estimate of the nodal error |y(rj,n) — Yo ee after 


<N,M 
extrapolation of YN", we decompose the solution Y"’“? on the mesh & ’ 
into the sum 
YNA@ _ RNAO 1 GN.Ae 


RN.A@ GN, A@ 


where the smooth component and the layer component are, 
respectively, the solutions of the following discrete problems: 
LNM RNAO — f in XNM RNAP —] on GNM 
eas . (24) 
LNMSN,A8 —Q in XVM SNAP m on GNM 
€ € ? ¢. d 


Y2N.46/2 into the smooth component 


—<2N,2M 
on the fine mesh & as 


Similarly, we decompose the solution 


R2N.48/2 and the layer component S$2-49/2 


Y2N,A0/2 _ p2N,A0/2 1 G2N,A0/2, 


Then, the errors can be written in the following form: 


YN Ae _ y= [ee = ‘ dhs Ga = m). 


4.3 Extrapolation of RN’? 


Let 


O71 


1 1 
mi(r, 0) = CH) and na(r, 0) = 2 Or2(r, 0)’ (r, 0) E X. 


Lemma 8. Assume that ¢ < N~!. Then, the local truncation error related 


to the smooth component fulfills the following condition: 
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Lae ee _ ‘) (73, On41) 
= hym (5, On41) + AOna(rj, On41) + o( op ae), forl<j<N-1. 


Proof. By utilizing the bounds on the derivatives of | as provided in Theorem 
2 and ¢ < N7! < GH, we can readily derive the following outcome from 


Taylor’s expansion: 


139 (16-1) 0un 


€ folal| foal h; 071 
= Al jpeg (1 On41) = hj 3 (w2,On41) + Waly) Fog (PI Am) 
he a] A6 8? Ab Bl, x 


= gy Y1(7s) ag (Ws, On) _ > age "i naa) = Sr gga "i? 5 


for some wy € (17,7 j41), W2, Ws E 500) and 6 E (An, On+i)- 


Now, following the classical approach of Keller [12], we define the functions 
Q;, where i = 1,2, as the solutions of the following IBVPs: 


LQ; =, in X, 
Qi(r,0) =0, reé Ay, (25) 
Q;(0,0) = Qi (1, 8) = 0, 0€ (0,71, 1= 1,2. 


Next, decompose Q; as Q; = 0;+®;, 7 = 1,2, where the smooth component 


OQ; and the layer component ©; satisfy the following IBVPs: 
LO; = Mh, L-®; =0 in X, 
O;(r,0) = ;(r,0) =0, re A,, 
0;(0,6) = (0,8) =0, O;(1,0) = —®;(1,6) 6€ (0,7), i=1,2. 
(26) 


Theorem 4. For all nonnegative integers p,q satisfying 0 < p+q < 2, the 
smooth components ©;,i = 1,2, defined in (26) satisfy the following bounds: 
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Proof. firstly, following the approach of [27] and using the bounds on the 
derivatives of | given in Theorem 2, we obtain that the solution Q; of the 
IBVP (25) satisfies the estimate 


oPtaQ; 
OrP 004 


(r, a) < o(ite exp(-a(l-n)/e), (r,0) EX, for 0< ptq < 2. 
(27) 


This shows that away from the side r = 1, the derivatives of Q; are bounded 


independent of ¢; see, for example, [33]. 


Remark 2. It is to be noted that the bounds obtained in Theorem 2 on 


P+] 

sans for 0<p+q< 4, are sufficient to derive the estimate (27) for the 
r 

solution Q; of the IBVP (25). 


Lemma 9. Assume that ¢ < N~!. Then, 


(ee = ' (75, On41) 
= h,Ox(r;, On41) + A003(7;, On+1) + O(n + a), for 1 < j < N-1. 


Proof. Assume that 1 < 7 < N —1. Then, applying Lemma 8 and (26), we 
get 


(x41) (rj, On+1 = hj L-O1(r;, On41)+AOL-Oa(r;, On4+1)+O (Hae) . 
(28) 


Now, following the Taylor’s expansion, it can be deduced that for 7 = 1, 2, 


é 830; 
(ae- 22") Oi(Ti, On41)] < Fa + hj+1) “Bre . 
hj 070; AO || 070; 
+ Burin) Or? || 7 2 | 00? a 


In accordance with Theorem 4 and the inequality h; < H for all j, where H = 
2</N, it follows that 


hj @ = i) Oi 78.5,) FAP @ = oe) 02(73,6n41) 


<c(H’+HA0+A0) <o(H’ +00"), (29) 


Hence, from the combination of (28) and (29), we get 
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Let [R*4"—1-nj0-000s| (aici ie o(H+a0), 1<j<N-1. 
(30) 


Provide the definitions for the following pair of discrete functions: 


UO (05; On) = C4 (w2+267) (l+rjJ4E(r;, (On), for 0 < j < Nandn > 0, 
where 


[RN:4? —1— hjO, — ACs] (rj,O,), forO<j<N-1, 

£(75, On) = 
0, for 7 =0,N. 

(31) 

Now, for 1 < j < N —1, using the inequality H < 2N~', where H = 2</N, 


we have 
Loe" ( (02406) (1+ 15)] > a(N-2+A6") > (H+ 00?) /4 (32) 


Therefore, by utilizing (30)and (32), one can select the constant C4 to be 
suitably large, ensuring that LN“ v+ > 0 within NY”. Additionally, in 


accordance with (24), (26), and (31), it is established that ut > 0 over 
GNM 
ae 


the domain Consequently, by employing Lemma 5 of the discrete 


minimum principle across the region Nl M we achieve 


| [RY 491-150, - 66, (rj, On41) 


ra o(w?+20"), forl<j<N-1, 


and hereby, the desired result follows. 


We deduce the estimate for the smooth part of the error after extrapola- 


tion, in the following lemma. 


Lemma 10. Assume that ¢ < N7~!. Then, the error after extrapolation 


RN:4?O 


associated to the smooth component satisfies 


U(75,On41) — Rovip (77s On41)| < o(w? a a6), for 1<j<N-1. 


Proof. Since the mesh widths of ail are half of those of ae applying 
—2N,2M 


Lemma 9 on the fine mesh ®, we have 
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(BAG, On41 — ') (r;, On+1) 


(FH /2)1 (175, On41) + (A0/2)02(15, On41) + o(w a a), 
for 1<7#< NP, 
at ee (33) 
(h/2)O1(175, On41) + (A/2)62(75, On41) + o(n-? - ae), 
for N/2+1<j<N-1, 


where H = 2¢/N and h = phv;_,. Therefore, according to the extrapolation 
formula (23), from Lemma 9 and (33), it immediately follows that 


U(rj, Anti) — Ree (Oy An+1) 


Ndi [2B AG, Baya) — BNA%(ry, an)) 
=——2 Cae — ') (Gan On+41) + (ae = ' (Ts On+1) 


=o(n? +00"), for 1<j<N-1. 


4.4 Extrapolation of SN’*? 


Lemma 11. The error after extrapolation associated to the layer component 


SNAt satisfies 


<ON~*, for 1<7<.N/2. 


ry On41) — Aaa Ce On+1) 


Proof. Let 1 < j < N/2. Now, following the argument resulting to [34] over 


yy and using (20), we can show that 


N ah —1 
<¢ |] (1+%) <CN~?. 


k=j+l 


s*4(5 On) 


Next, from (19) and Theorem 2 we obtain |m(rj, On+41| < CN~?. Hence, we 


have 
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| ee am m) (re; On41) < EN. 


SON Au? 2 m) (rj, On41)| < CN~? and hereby, the required 


Similarly, | ( 


extrapolated error bound is obtained from the extrapolation formula (23). 


Qe 02m 10?m 
‘ Let ai(r, 0) = @ ir) aa lr 6), and let g(r, 6) = 2 092 0), (r,O) € 


Lemma 12. For N/2+1< 7 < N-—1, the local truncation error associated 


to the layer component satisfies 


rat Ge = m) (ee Baa) = (w- in(1/e)) 01(7;, On41) + A002(7;,9n+1) 


+O G exp(—a(1 — r341)/e)N7? In? (1/e) 


+exp(—a(1—1;)/e) A #) : 


Proof. Let N/2+1 < j < N-—1, Then, from the Taylor’s expansion and using 


Theorem 2, for some w) € (1;,7541),W2,W3 € (rj-1,7;), and Oe (An; On+1); 


we have 
LNM Ga _ m) Cath 
= ae Peon a aL + Pyar) (rj Oat) 
= Pw (r5) (25. On) + A or Bn) - Oe oy, p) 


De fo 02m Ab 02m 
= (w “in(d/e) ) Wale) Sep (roma) + > gz (tir tnt) 


+O G exp(—a(1 — r341)/e)N~? In?(1/e) + exp(—a(1 — r;)/e) A #). 


Let the functions Q;,2 = 1,2, represent the solutions of the following 
IBVPs: 


Iran. J. Numer. Anal. Optim., Vol. 14, No. 1, 2024, pp 219-264 


Sah and Gowrisankar 242 


L-Q@ =o in Q= (1 —€é, 1) x (0, 7], 
Q (7,0) =0, re fl—e,]], (34) 
Q(1—¢,0) =Q,(1,6)=0, 6€(0,T7], i=1,2. 

The bounds on the derivatives of Q;,i = 1,2, are obtained by using Lemmas 


13-14. We provide a detailed proof for the function Q; and show that the 


proof for Q2 follows in a similar manner. 


Lemma 13. The bounds for the temporal derivatives of Q; are as follows: 


7Q1 
eq 0)| < Cexp(—a(1l—r)/e), (r,0)€Q, for 0<q <3. 


Proof. From (34), we have Q; =0 on Teg =2\Q. Also, from Theorem 2, 


we obtain 
|L-Qi(r, 9)| = |o1| < Ce~* exp(—a(1—r)/e), (r,0) €Q. 
Now, consider the barrier function 
w(r, 0) = Cexp(—a(1—1)/e), (7,6) €Q, (35) 


where C is sufficiently large. It follows that Lew(r,0) > |L-eQi(r, @)| in 2 
and w > Q, on Tq. Since Lz satisfies the minimum principle Lemma 1 on 
Q, there we have Q; < Cexp(—a(1—r)/e), in Q. 


Next, we shall obtain the bound on On the sides r = 1 — € and 


a O 
r = 1 of Q, we have Q; = 0 and so = 0. On the side 0 = 0, we have 
O o? 
Q; = 0, and so Qi = @ = 0. Hence, from (34), we obtain 
Or Or? 
OQ) 2€ Oem 
59 (79) = G eur) aa (9), re (1 —e,1]. (36) 
2m 
Since on the side 6 = 0, we have m = 0, which implies that 5,2 = 0, and 
ir 
hence from (36), we obtain or, 0) =0, re [1—e,1]. Therefore, it is clear 
that 
OQ 


59 0) = 0, on Ta. 


Now, differentiating the equation given in (34) with respect to 0, we get 
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Qi PQ, 07Q, OQi _ 2 am 
ae abare + M1) aga, + M2 Ge = TM) aga 2 


in Q, (37) 


and employing Theorem 2 in (37) yields that 


1.(S) (r, 0) < Cel exp(—a(1—r)/e), (7,0) € 2. 


Therefore, choosing the barrier function 7(r,@) given in (20), we obtain that 


LewW(r,0) > |Le (FF) (r,8)| in Q and w > on Tq. 
Thus, applying the minimum principle on 2, we have 
0 
Qi < Cexp(—a(l—r)/e) in 2. 
00 
ie PQ 
Similarly, we shall obtain the bound on Fee On the side r= 1—e andr = 
0 of Q Hence, 
fag 
str, 8) < Cexp(—a(1—r)/e), (r,@) € ©. 
Finally, by employing a similar approach as for the second-order derivative 
2 3 
ae one can derive the bound for a thereby completing the proof. 


Lemma 14. The following mixed derivative of Q;, bound on the derivative 
of a and spatial derivative of Q, satisfies the following bound: 


(A). a" (r, a < Ce~exp(—a(1—r)/e), (7,0) € 


Or? - 


(B). “ (SF) fa) < Ce? exp(—a(l—r)/e), (r,0) € Q for p = 1,2. 


aPQu = 
(C). Aye (r, 8) < Ce exp(—a(1—r)/e), (7,0) €Q forl<p<3. 


Proof. For (A) from Lemma 13, we have 


Of PQ FQ PQ 
<5 (rage) PMU) pagge PUA) pa ee ES ABE) 
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OQ, 2€ o*m 
e998 gt) poner 


where Ji(r, 0) = 


2, we have 


From Lemma 13 and Theorem 


|Ji(r,0)| < Ce~* exp(—a(1—r)/e), (r,0) € Q. (39) 


3 

Let 6 € [0,T]. be fixed. To obtain the desired bound on a we apply the 
r 

argument presented by Kellogg and Tsan [13] on the line segment {(r,@),r € 


[1 — ¢,1]}. This requires to use (38), inequality (39), and the previously 
PQ 
06? ° 


established bound on 


For (B) from Lemma 13, rewriting (37), we get the following form: 


0 Qi 0?Q, 5 OQ, = 
-<5: (FS) + Wir) ae + Wo(r) 5 =Jn EQ, (40) 
2 3 
Ja(r, 6) = ge Gs + 2 WA (r) oe From Lemmas 13 and 14 (A) and 


06? a 000r2" 
Theorem 2, then we get 


OP Jo 
OrP 


(r, a) < Cel exp(—a(1—r)/e), (7,0) € Q. (41) 


Now, let us consider the fixed time @ € [0,7]. Applying the argument pre- 
sented by Kellogg and Tsan [13] to the line segment {(r,0),r € [1 — ¢, l]}, 


using (40) inequality (40) and the bound on ae we can obtain the desired 
2 
Similarly, differentiating (40) with respect to r and use of 
i; 
6) o? 
the inequality obtained bound on a 
1 


Or200° 


bound for 


directly implies the 


desired bound on 


For (C) from Lemma 12, rewriting (34) we get the following form: 


Od (AQ) | OQi 7 
5 ( a ) T Wi (r) ar + WoQ, = J3 €D, (42) 
2 
where J3(r,@) = OQ: | 2. y dg “From Lemmas 13 and 14 (B), we have 


00 ' a’ Or2 
OP Js 
OrP 


(r; a) < Ce exp(—a(1—r)/e), (r,0)€Q for O<p<2. (43) 
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Let us now fix @ € [0,7] and proceed to apply the argument presented by 
Kellogg and Tsan [13] to the line segment {(r,6),r € [1 — ¢,1]}. By utilizing 


(42), inequality (43), and the given bound on Q:, we can derive the desired 


bound for a , 


The bound on = 


can be deduced in a similar manner as in the previous 


201 
Or 
OT aQ 


establish the desired bound for on Lastly, to obtain the bound on = % 
r ir 


one can differentiate with respect to r and employ (42) in a similar fashion. 


cases. By employing (43), along with the bounds on Q; and , we can 


Subsequently, we consolidate the preceding lemmas into the following the- 


orem. 
Theorem 5. For all nonnegative integers p,q, satisfying 0 < p+q < 3, the 
solution Q;,i = 1,2, of (34) satisfies the following bounds: 


artig, 
OrP 4 


(nr, a) < Ce? exp(—a(1—r)/e), (r,0) €Q. 
Lemma 15. For N/2+1< j< N—1, we obtain 

Gig = m) (rj, On+1) = (w- in(1/e) ) Qu(r, 8n+1) + MOQ2(r;,9n41) 
+ o(x-?nt(a/e) - a6). 


Proof. Let N/2+1<j<N-1. 
Then, Lemma 12 and (34) imply that 


133 (s84 mn) rn 
— (w> in(1/e)) LQilr;, On41) + AOL -Qa(r;, On41) 


+O (= exp(—a(1 — r)/e)N~? In?(1/e) + exp(—a(1 — r;)/e) A *) 
(44) 


Again, from the Taylor’s expansion and Theorem 5 it follows that for / = 1, 2, 


@ -_ 7 Q;(r;, On+1) 
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= o(e exp(—a(1 — 1r541)/e)N~' n(1/e) + exp(—a(1 — r;)/e) A 6), 
and so 


| (w- in(1/e) @ - 188) Ours, On4+1) + AO @ - LE) Oats, On+1) 


< cle exp(—a(1 — rj41)/e)N~?In-?(1/e) + exp(—a(1 — r;)/e) A 6? 
+ G exp(—a(1 — rj41)/e) + exp(—a(1 — rie) (w- in(1/e)) A ” 


< Ce! exp(—a(1 — rj41)/e) (w In?(1/e) + 6) : 


(45) 
Therefore, combining (44) and (45), we have 
alg sve —m- (w- In(1/e)) - 005| (Ti, On41) 
<Cs Ge exp(—a(1 — 7541)/e) (w In?(1/e) + As") ), 
for N/2+1<j<N-1, (46) 


for some constant C5. Now, consider the barrier function 


N —1 
hy 
xP =O [wa 4+r;)+ (-?m2a4/e) - Ae?) Ry II (1 + =~) | 
k=1 


for N/2<j<N. 


For N/2+1<j< N—1, employing (18), we get 


N -1 
h 
LIM Gt" > CoCee* (1-2 m2(4/e) + a6) I] (: ey 


k=j+1 
Also, for 7 > N/2, from (19), we get 
N ah —1 

exp(-all —ryai)/e) Sexp(anje) T] (2+ 4) 

k=j+1 . 

N ah -1 

<exp(r) [] (1 + a) (48) 
k=j+1 : 
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Hence, for N/2+1<j< N-—1, it is easy to show from (46) — (48) that 


PO a a aaa, eee |svee —m— (w- In(1/e) ) Q — 60s| (73, On+1)]; 
(49) 


provided Cg > exp(A1)Cs/C2. Again, from the proof of Lemma 11, we get 


| (snes = m) (l= fea) 


< C7N~?, 


for some constant C7. Hence, by (34), 


VeteN eS | ee =—m— (w- In(1/e)) = 00s| (i-efe34)\, 
(50) 


provided Cg > C7. Also, recalling (24) and (34), we have 


ye SS 0= | isvne —-m- (w- In(1/e)) - 605| (1, Prada). (1) 


Therefore, it follows from (49)—(51) that one can choose 
Co = max{exp(A1)C5/Co, Cr} 


so that x7 is a barrier function for 


+ oe —-m- (w- In(1/e)) _ 605| (73,On)- 


Thus, by the discrete minimum principle over aaa 1) ([1 —€, 1] x [0, T Ds 


we have 


| ane —m-— (w= in(1/e)) Q = 60s| (73; On41) 


41 
SXF 


= o(w In?(1/e) + a0), 


for N/2+1<¥7< N-—1, and this completes the proof. 


Lemma 16. The error after extrapolation associated to the layer component 


SN:4?9 satisfies 
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ry Patt) = Cas Ona) < o(w? In?(1/e) + a), 


for N/2+1<j<N-1. 
Proof. Let N/2+1<j7< N-—1. from Lemma 15, we have 


(s*4"—m) (75, Ont1) = N~*Q1 (175, On41)+AOQa (15, On41)+O( N~?+A6? J. 
(52) 
: <2N,2M re ' 
Next, since the fine mesh XN, has the same transition point 1—e as that 


of ee Lemma 15 implies that 
(ener = m) (5, On41) 
= (2N)Qi(r;, On41) + AO/2Q2(15, On41) + O (w + 6) ° (53) 


Henceforth, eliminating O(N~!) and O(A@) terms from (52) and (53), the 


required extrapolated error bound is obtained. 


4.5 Convergence result of the solution pial el 


The main result of this paper is presented in the following theorem. 


Theorem 6. [Error after extrapolation] Assume that « < N~!. Let y be the 


solution of the continuous problem (1) — (2) and let ye be the solution 


obtained via the Richardson extrapolation technique by solving the discrete 
~<2N,2M 


problem (16) on two nested meshes RM and X. Then, the error 


associated with the solution arg ® at time level 0,, satisfies 


N,AO 
yrs, On) > Voges (rj, On) 


= o(w-?na/e) + a6), for 1<j<N-1. 
(54) 


Proof. For each (rj, On) € cae we have 


ulr,9n) — ¥N2? = (ur), 6,) - a) + Ge 6,) - sit) , 
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Therefore, the result (54) follows immediately after combining Lemma 10 for 


the smooth component and Lemmas 11 and (16 for the layer component. 


5 Numerical results 


In this section, we showcase the numerical findings acquired through the 
utilization of the Richardson extrapolation technique, aimed at corroborating 
the theoretical outcomes asserted in the preceding section. To accomplish this 
objective, we conduct extensive numerical experiments on a modified graded 
mesh ae For all the test examples, we carefully select the transition 
parameter p = 0.9, along with a = 1 anda time step size of AO = 1/N. These 


parameters are consistently applied throughout the numerical investigations. 


Example 1. Let us examine the parabolic IBVP presented below: 


dy & a 
ag epee + (+r) = Fr, 8), 707.8) € (0,1) x (0,4), 
y(r,0) =yo(r), O<r<l, (55) 


y(0,0) =0, y(1,0)=0, 0<0<1. 


Example 2. Let us examine the parabolic IBVP presented below: 


Oy O*y dy 
36 fae t Bp TL), (r,) € (0,1) x (0, 1], 
y(r, 0) = yo(r), O0<r<l, (56) 


y(0,0) =0, y(1,@) =0, 0<¢<1, 


where the initial data yo(r) and the source function f(r,@) have been 
chosen to fit the exact solution for the IBVP (55) and (56), which is given 


by the expression 


y(r, 0) = exp(—6) (mi + 72r — exp (- : = “)) 


where the constants are defined as m1 = exp (-+) and m2 = 1 — exp (—-2) i 


Once the exact solution is known for each ¢, the maximum point-wise error 


iN before and after extrapolation is calculated using 
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( onjeRY M y(r5, On) aes (rj, a) 
75 ,On JER.” 
and 
N,A@ 
( sieRY M y(T5, On) — vetin (rj, On) ’ 
75 ,On)ERe? 


respectively, where y(r;,6n), bad ar 6,) and Yo Caen) respectively, 


denote the exact solution, the upwind numerical solution and the extrapo- 
lated solution obtained on the mesh nM with N mesh-points in the spatial 
direction and M mesh-points in the @-direction such that A@ = T’/M is the 
uniform time step. In addition, the corresponding order of convergence is 


determined by 


=N,A@ 
py .49 =]o 4E 
€ = 1082 2N,A0/2 }° 
—E 


Now, for each N and AQ, we define Fa A? — max, a 9 as the e-uniform 


maximum point-wise error and the corresponding local ¢-uniform order of 


convergence is defined by 


N,A@ 
PN“? = log, | Saas 
= 1082 | Gon Ao2 }° 
>N,AO 


Tables 1 and 3 present the calculated maximum point-wise errors Ez 


, as 


well as the corresponding order of convergence P’:“? for various values of € 
and N before and after extrapolation, respectively, for Example 1 and (55). 
Similarly, Tables 2 and 4 showcase the corresponding results for Example 2 


and (56). The results presented in Tables 1-4 demonstrate a consistent de- 


crease in the computed ¢-uniform errors EN“? for Examples 1 and 2 as the 


number of grid points, N, increases. Besides, the comparison of numerical 
results obtained by the proposed scheme and results in [18] are tabulated in 
Tables 5 and 6 for Example 1. From these tables, one can conclude that the 
proposed scheme gives better results than the scheme considered in [18]. Fig- 
ures 2—5 exhibit the solution profile for Examples 1 and 2 for various values of 
eand N. Figures 2 and 3 provide a solution profile before Richardson extrap- 
olation whereas Figures 4 and 5 provide a solution profile after Richardson ex- 
trapolation. This observation indicates that the implicit upwind scheme (16) 
exhibits ¢-uniform convergence before and after extrapolation. Additionally, 


Figure 6 depicts the maximum point-wise errors for Examples 1 and 2. The 
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plots clearly illustrate that Richardson extrapolation enhances the order of 
convergence of the upwind scheme from O(N~1In(1/e)) to O(N~? In? (1/e)) 
for A@ = 1/N. This finding validates the theoretical bounds established in 
(16) and (33). 


6 Conclusion 


In this article, our focus revolved around devising an efficient numerical ap- 
proach for tackling one-dimensional singularly perturbed parabolic problems. 
Specifically, we leveraged the Richardson extrapolation technique to enhance 
the performance of the upwind scheme in addressing the one-dimensional sin- 
gularly perturbed parabolic convection-diffusion IBVP presented in (1) — (2). 
The methodology began by discretizing the domain using a modified graded 
mesh. This mesh structure was defined by N x M points. Subsequently, 
we applied the classical implicit upwind finite difference scheme to solve the 
IBVP effectively. To validate the efficacy of our approach, we furnished rigor- 
ous proof detailing the ¢-uniform error estimation. This proof substantiated 
the convergence of the upwind scheme in an ¢-uniform manner, displaying 
nearly first-order accuracy. Furthermore, by maintaining the transition pa- 
rameter at a constant value of p = 0.9 within the modified graded mesh 
framework, we extended our analysis to the fine mesh configuration, which is 
characterized by 2N x 2M mesh-points. Finally, we amalgamated the out- 
comes obtained from both coarse and fine mesh computations. This amal- 
gamation served as a foundation for the implementation of the Richardson 
extrapolation technique. Both theoretical and computational examinations 
affirmed that this extrapolation technique propels the convergence of the 
basic upwinding scheme from almost first-order accuracy to a nearly second- 


order accuracy, all within the confines of the discrete supremum norm. 


Iran. J. Numer. Anal. Optim., Vol. 14, No. 1, 2024, pp 219-264 


Sah and Gowrisankar 252 


(c) (d) 


Figure 2: Solution profile for Example 1 using simple upwind scheme before Richardson 
extrapolation technique with different value of ¢ and N. (a) Solution of N = 128,e= 107? 
(b) Solution of N = 128,e=10~* (c) Solution of N = 128,e= 107° (d) Solution of 
N = 128,e= 107-8. 
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(c) (d) 


Figure 3: Solution profile for Example 2 using simple upwind scheme before Richardson 
extrapolation technique with different value of ¢ and N. (a) Solution of N = 128,e= 107? 
(b) Solution of N = 128,e=10~* (c) Solution of N = 128,e= 107° (d) Solution of 
N = 128,e= 107-8. 
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Figure 4: Solution profile for Example 1 using simple upwind scheme and after 
Richardson extrapolation technique with different value of « and N. (a) Solution of 
N = 128,e= 10~? (b) Solution of N = 128,e= 10-4 (c) Solution of N = 128,e= 10-6 
(d) Solution of N = 128,e= 107-8. 
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(c) (d) 


Figure 5: Solution profile for Example 2 using simple upwind scheme and after 
Richardson extrapolation technique with different value of « and N. (a) Solution of 
N = 128,e= 10-2 (b) Solution of N = 128,e= 10-4 (c) Solution of N = 128,e= 10-6 
(d) Solution of N = 128,e= 10-8. 
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Table 1: Maximum point-wise errors and the corresponding order of convergence before 


extrapolation for Example 1 on Modified graded mesh 


Number of Intervals N/A@ 


€ 128/75 256/35 512/4 1024/35 2048/ hy 4096 / 355 


07? 1.3707e — 02 7.5039e — 03 3.9359e — 03 2.0188e — 03 1.023le — 03 5.1505e — 04 


0.8692 0.9310 0.9632 0.9805 0.9902 


Or 3.2822e — 02 1.7760e — 02 9.1788e — 03 4.7257e — 03 2.3904e — 03 1.2045e — 03 


0.8860 0.9523 0.9578 0.9832 0.9889 


0-8 4.8001le — 02 2.6420e — 02 1.3866e — 02 7.1029e — 03 3.5934e — 03 1.8063e — 03 


0.8614 0.9302 0.9650 0.9831 0.9923 


o-8 6.2334e — 02 3.4617e — 02 1.8304e — 02 9.4150e — 03 4.7736e — 03 2.4034e — 03 


0.8485 0.9193 0.9591 0.9799 0.9900 


10-1 7.4394e — 02 4.2641e — 02 2.2683¢e — 02 1.1706e — 02 5.9463e — 03 2.9968e — 03 


0.8030 0.9106 0.9544 0.9772 0.9886 


EN.Aé 7.4394e — 02 4.264le—02  2.2683e — 02 1.1706e —02 5.9463e-—03 —-2.9968e — 03 


pN.Aée 0.8030 0.9106 0.9544 0.9772 0.9886 


Table 2: Maximum point-wise errors and the corresponding order of convergence before 


extrapolation for Example 2 on Modified graded mesh 


Number of Intervals N/A@ 


128/44 256/4 512/4 1024/4 2048/34, 4096/45 


™ 


0-? 1.3635e — 02 7.4761e — 03 3.9204e — 03 2.0104e — 03 1.0185e — 03 5.1265e — 04 


0.8669 0.9313 0.9636 0.9811 0.9903 


o-4 3.2821e — 02 1.776le — 02 9.1790e — 03 4.7258 — 03 2.3905e — 03 1.2045e — 03 


0.8859 0.9523 0.9578 0.9832 0.9889 


0-8 4.7999e — 02 2.6420e — 02 1.3865e — 02 7.1029e — 03 3.5934e — 03 1.8063e — 03 


0.8613 0.9302 0.9650 0.9831 0.9923 


o-8 6.2332e — 02 3.4617e — 02 1.8304e — 02 9.4150e — 03 4.7736e — 03 2.4034e — 03 


0.8485 0.9193 0.9591 0.9799 0.9900 


10-1 7.4392e — 02 4.2641e — 02 2.2683¢e — 02 1.1706e — 02 5.9463e — 03 2.9968e — 03 


0.8029 0.9106 0.9544 0.9772 0.9886 


ENA 7.4392e — 02 4.2641e — 02 2.2683¢e — 02 1.1706e — 02 5.9463e — 03 2.9968e — 03 


PN,A8 0.8029 0.9106 0.9544 0.9772 0.9886 
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Table 3: Maximum point-wise errors and the corresponding order of convergence after 


extrapolation for Example 1 on Modified graded mesh 


Number of Intervals N/A@ 


m 


128/%5 256/35 512/ a 1024/35 2048/ sty 4096/ 355 


10-? 1.8300e — 03 6.7491le — 04 2.1138e — 04 6.1563e — 05 1.6744e — 05 4.3943e — 06 


1.4390 1.6748 1.7797 1.8784 1.9299 


10-4 1.3132e — 03 4.8091le — 04 2.4003e — 04 1.6554e — 04 1.2649e — 04 8.9042e — 05 


1.4492 1.0025 1.5360 1.3881 1.5065 


10-6 2.5422e — 03 7.2236e — 04 1.9374e — 04 5.1776e — 05 1.4992e — 05 5.6506e — 06 


1.8153 1.8986 1.9038 1.7881 1.4077 


10-8 4.361le — 03 1.2692e — 03 3.3970e — 04 8.8050e — 05 2.2424e — 05 5.6808e — 06 


1.7808 1.9015 1.9479 1.9733 1.9809 


19710 6.5786e — 03 1.9550e — 03 5.2916e — 04 1.3762 — 04 3.5092e — 05 8.8610e — 06 


1.7506 1.8854 1.9430 1.9715 1.9856 


BNAG 6.5786e — 03 1.9550e — 03 5.2916e — 0¢ 1.3762e — 04 3.5092e — 05 8.8610e — 06 


pN.A@ 1.7506 1.8854 1.9430 1.9715 1.9856 


Table 4: Maximum point-wise errors and the corresponding order of convergence after 


extrapolation for Example 2 on Modified graded mesh 


Number of Intervals N/A@ 


€ 28/45 256/35 512/ a 1024/5 2048/ sty 4096 / 335 


10-? 1.8144e — 03 6.6235e — 04 2.0693e — 04 6.0125e — 05 1.6335¢e — 05 4.2857e — 06 


4538 1.6784 7831 1.8800 -9304 


10-4 1.3142e — 03 4.8110e — 04 2.4012e — 04 1.6561e — 04 1.2653e — 04 8.9056e — 05 


4498 1.0026 5360 1.3882 5067 


10~ 2.5430e — 03 7.2242e — 04 1.9375e — 04 5.1777e — 05 1.4991e — 05 5.6505e — 06 


a 


8156 1.8987 9038 1.7882 A077 


10-8 4.3619e — 03 1.2692e — 03 3.3971le — 04 8.8050e — 05 2.2424e — 05 5.6808e — 06 


-7810 1.9016 9479 1.9733 -9809 


197° 6.5793e — 03 1.9550e — 03 5.2917e — 04 1.3762e — 04 3.5092e — 05 8.8617e — 06 


-T507 1.8854 9430 1.9715 -9855 


BNAG 6.5793e — 03 1.9550e — 03 5.2917e — 04 1.3762e — 04 3.5092e — 05 8.8617e — 06 


pN.A@ 1.7507 1.8854 1.9430 1.9715 1.9855 
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Table 5: Comparison of maximum point-wise errors and the corresponding order of 


convergence before extrapolation for Example 1 


Number of Intervals N/A@ 


é 128/4 256/35 512/% 1024/3, 2048/7455 4096/35 

10-2 1.3707e —02 7.5039e —03  3.9359e —03  2.0188e 03  1.023le—03 —5.1505e — 04 
0.8692 0.9310 0.9632 0.9805 0.9902 

10-4 3.2822e —02 1.7760e-02 9.1788e-03 4.7257e-03 2.3904e-03  1.2045e — 03 
0.8860 0.9523 0.9578 0.9832 0.9889 

10-6 4.800le —02 2.6420e-—02 1.3866e-02 7.1029°-03 3.5934e-03 —1.8063e — 03 
0.8614 0.9302 0.9650 0.9831 0.9923 

EN.Ae 4.800le —02 2.6420e 02 1.3866e-02 7.1029°-03 3.5934e-03 —1.8063e — 03 
PN.Aa 0.8614 0.9302 0.9650 0.9831 0.9923 

Results in [18] 
Number of Intervals NV 
€ 32 64 128 256 512 1024 

10-? 4.3463e —02 2.9874e-02 1.876le—02 1.1330e—02 6.5980e—03  3.7392e — 03 
0.5408 0.6711 0.7276 0.7799 0.8192 

10-4 4.2785e —02 2.9447e-02 1.8513e-02 1.1186e—02 6.5158e—03  3.6933e — 03 
0.5389 0.6695 0.7269 0.7796 0.8190 

10-® 4.2779e —02 2.9443e-02 1.851le—02 1.1184e-02 6.5150e-03  3.6928e — 03 
0.5389 0.6695 0.7269 0.7796 0.8190 

ENA 4.2779e —02 2.9443e 02 1.851le—02 1.1184e—02 6.5150e—03 3.6928 — 03 
pN,48 0.5389 0.6695 0.7269 0.7796 0.8190 
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Table 6: Comparison of maximum point-wise errors and the corresponding order of 


convergence after extrapolation for Example 1 


Number of Intervals N/A@ 


é 128/45 256/35 512/45 1024/5 2048/at5 4096/3355 

10-2 1.8300e —03  6.7491e—04 2.1138e—04 6.1563e-—05 1.6744e—05 4.3943e — 06 
.4390 1.6748 1.7797 1.8784 1.9299 

10-4 1.3132e —03  4.8091le—04 2.4003e—04 1.6554e-04 1.2649 — 04 8.9042e — 05 
.4492 1.0025 1.5360 1.3881 1.5065 

10-8 2.5422e —03 7.2236e 04 1.9374e-04 5.1776e—05 1.4992e-05 5.6506e — 06 
8153 1.8986 1.9038 1.7881 1.4077 

ENA? 2.5422e —03 7.2236e-04 1.9374e-04 5.1776e-05 1.4992e -—05 5.6506 — 06 
pN.Ae 1.8153 1.8986 1.9038 1.7881 1.4077 


Results in [18] 


Number of Intervals NV 


€ 32 64 128 256 512 1024 

10-2 6.7626e —03 2.9552e 03 1.2260e-03 4.4704e—03  1.5112e-—04 4.8606¢e — 05 
1.1943 1.2693 A555 1.5646 1.6365 

10-4 6.7299 —03 2.930803 1.2146e—03 4.4264e—04 1.4950e—04 4.8044e — 05 
1.1993 1.2708 4563 1.5660 1.6377 

10-6 6.7297e —03 2.9306e 03 1.2145e-03 4.4260e-04 1.4949e 04 —4.8039e — 05 
1.1993 1.2708 4563 1.5660 1.6378 

EN.A@ 6.7297e —03 2.9306e 03 1.2145e-03 4.4260e-04 1.4949e —04 —4.8039e — 05 
pN.Ae 1.1993 1.2708 4563 1.5660 1.6378 
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Figure 6: Before extrapolation and after extrapolation Loglog plot of the maximum 
point-wise errors for Example 1 and Example 2. (a) Before extrapolation for Example 1 


(b) Before extrapolation for Example 2 (c) After extrapolation for Example 1 (d) After 


extrapolation for Example 2. 
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